A Unified Spectroscopic and Photometric Model to Infer Surface Inhomogeneity: Application to Luhman 16B

M dwarf stars comprise the most numerous spectral class in the Milky Way. Due to M dwarfs' small size, they are attractive targets for studies attempting to find habitable terrestrial planets because orbiting planets create correspondingly larger and more easily detected signatures with both the transit and radial velocity (RV) methods. Small rocky planets may be over 3 times more common around M dwarfs than FGK stars (Mulders et al. 2015). The overabundance was later confirmed in comparing RV and transiting surveys (Tuomi et al. 2019; Sabotta et al. 2021; Zink et al. 2020).

Despite the prevalence of planets around M stars, obstacles exist for declaring M dwarf stellar systems ideal for planetary habitability. M dwarfs' low temperatures (compared to G dwarf stars like our Sun), mean planets located at distances where liquid water could exist on a planet's surface, approximately 10% the distance between Earth and the Sun, would experience extensive stellar activity including flares (Kasting et al. 1993; Kopparapu et al. 2013). Fast-rotating M dwarfs, and particularly cooler, later and fully convective stars, are far more likely to flare than their hotter counterparts, which possess radiative zones (Günther et al. 2020). Furthermore, flares and even superflares have been observed in ultracool objects with spectral classes as cool as L5, demonstrating that intense magnetic activity continues beyond the M-class range (Paudel et al. 2018, 2020).

Flares and coronal mass ejections associated with star spots could potentially bathe nearby, likely tidally locked, planets in X-ray and ultraviolet (UV) radiation and charged particles, stripping the planets' nascent atmospheres. However, it remains uncertain whether this radiation will inhibit or catalyze the emergence and growth of burgeoning life (Scalo et al. 2007; Tarter et al. 2007; Zendejas et al. 2010). Laboratory experiments suggest that stellar flares could provide the UV light necessary to photochemically produce the prerequisite molecules for abiogenesis (Rimmer et al. 2018).

Understanding starspot coverage as a proxy for magnetic activity and stellar flares is necessary to reveal how these energetic events may influence the long-term habitability of planetary systems. Late M dwarfs have been observed to have star spots preferentially located within their polar latitudes which differs from the predominantly equatorial-forming Sun spots (Strassmeier & Rice 1998; Järvinen et al. 2008; Barnes et al. 2015, 2017). TESS observations have also provided preliminary evidence of white-light flares preferentially occurring at polar latitudes (between 55° and 81°) in fully convective M dwarf stars (Ilin et al. 2021). Future M dwarf imaging will be critical to confirming theories on starspot formation and migration, and will also be necessary to understand and quantify habitability in these star systems.

Lower in mass than M dwarfs, BDs are substellar objects below the hydrogen-burning minimum mass (M ≲ 70 M_J ) (Chabrier et al. 2000; Dupuy & Liu 2017). They gradually cool over billions of years from early post-formation temperatures consistent with late M dwarfs, through L, T, and Y substellar classes. Across this span, the spectral features transition quite dramatically in J − K color from red, cloudy L Dwarfs to relatively blue, cloud-free mid-T Dwarfs due to potentially two primary mechanisms related to the decreasing temperature: (1) Clouds composed of refractory materials such as iron (Fe), corundum (Al₂O₃), enstatite (MgSiO₃), and forsterite (Mg₂SiO₄) condense and precipitate out over the L/T transition leading to these materials seemingly being sequestered below the photosphere (2) carbon monoxide (CO) is reduced to methane (CH₄) resulting in J-band brightening and H and K suppression due to the ensuing CH₄ absorption (Tsuji et al. 1996; Jones & Tsuji 1997; Noll et al. 2000; Allard et al. 2001; Marley et al. 2002; Burrows et al. 2006; Cushing et al. 2008).

BDs at the L/T transition appear to demonstrate cloud patchiness, which may explain the source of variability. It was initially theorized that the condensation of the metal and silicate clouds along with robust tropospheric convection leads to gaps in cloud coverage (Ackerman & Marley 2001; Burgasser et al. 2002; Marley et al. 2010). More recently, it has been argued that a more likely scenario is that condensation results in patchy regions of varying cloud height and optical thickness (Radigan et al. 2012; Apai et al. 2013; Buenzli et al. 2014). Inhomogeneous atmospheres paired with the fact most BDs spin rapidly with average rotational speeds exceeding 20 km s⁻¹ for mid-L dwarfs, naturally leads to the conclusion that ultracool objects at the L/T transition will experience large photometric and spectroscopic variability (Reiners & Basri 2008).

Previous observations indicated that, although objects at the L/T transition were variable, they were perhaps no more variable than other BDs (Enoch et al. 2003; Buenzli et al. 2014; Metchev et al. 2015). However, Radigan et al. (2014) deduced that ${39}_{-14}^{+16} \%$ of L/T transition stars were highly variable (>2%) while also demonstrating strong variability to be rare for BDs outside this range. Eriksson et al. (2019) determined a similar value of ${40}_{-19}^{+32} \%$ for the proportion of strongly variable L/T transition BDs.

Future opportunities exist to shed more light on the exact nature of the L/T transition and BD atmospheric dynamics. ELTs and the JWST will provide high quality data that could potentially be used to map the surfaces of numerous nearby BDs. These maps very well may reveal the existence of atmospheric structure including zonal jets and vortices as well as the true source of L/T transition variability. With this information, we can further constrain atmospheric formation, structure, and dynamics.

Not only will we be able to map these ultracool objects' surface temperatures, we may also be able to map chemical tracers via Doppler imaging. This could allow us detect regions of chemical disequilibrium (Section 5.3) originating from vertical atmospheric flow. Studying such behavior across spectral classes can allow us to fully understand the development of thermal and hydrodynamic structures with ultracool objects.

Directly imaged gas giant exoplanets typically form at large distances from their host stars and are weakly illuminated, leading to targets similar to BDs in terms of observational opportunities. Like their more massive BD cousins, gas giant exoplanets possess clouds (Madhusudhan et al. 2011; Currie et al. 2011, 2014; Marley et al. 2012; Skemer et al. 2014) and chemical disequilibrium (Barman et al. 2011, 2015; Konopacky et al. 2013; Skemer et al. 2014; Currie et al. 2014; Lavie et al. 2017; Mollière et al. 2020; Wang et al. 2020).

High S/N and high spectral resolution data from ELTs will be essential to creating gas giant surface maps via Doppler imaging in the future. These maps may have the ability to determine the extent of cloud patchiness as well as the occurrence of planetary-scale structures such as zonal jet bands, vortices, and persistent storms. As is later discussed in Section 5.3, Doppler imaging may also be used to identify regions of convective activity through chemical disequilibrium tracers.

Doppler imagining is an important spectroscopic method for studying surface inhomogeneity. The history of Doppler imaging goes back almost 70 yr when Deutsch (1958) first used spectral line profile (LP) deviations to map the magnetic fields and abundance anomalies of Ap stars; this method was further formalized in Deutsch (1970) and Falk & Wehlau (1974). Ap stars continued to be targeted for surveys mapping their spots and chemical compositions due the class possessing strong surface inhomogeneities with century-long lifetimes (Khokhlova 1976; Goncharskii et al. 1977, 1982). Following this work, Vogt & Penrod (1983) conducted a groundbreaking study involving the Doppler imaging of the spotted star, HR1099.

These early Doppler imaging efforts suffered from the fact that stellar surface imaging is an ill-posed problem, the results were often degenerate with more than one solution for a given set of data (Piskunov et al. 1990). Solutions to this problem arose in applying maximum entropy methods (Vogt et al. 1987) and Tikhonov regularization (Piskunov et al. 1990) to infer the smoothest stellar surfaces.

Since these developments, studies have generated Doppler images of giants and sub-giants (Strassmeier et al. 1998, 1999; Donati 1999; Strassmeier 1999), pre-main-sequence stars (Hatzes 1995; Strassmeier & Rice 1998; Finociety et al. 2021), K Dwarfs (Collier Cameron & Unruh 1994; Barnes et al. 2001, 2004; Zaire et al. 2021), and even M dwarfs (Barnes & Collier Cameron 2001; Barnes et al. 2004, 2017). Notably, Crossfield et al. (2014) conducted the first Doppler imaging of a BD on the nearby binary system Luhman 16; these results were reprocessed and replicated in Luger et al. (2021).

Photometric and low-resolution spectroscopic (spectrophotometric) campaigns have led to improved understanding of BD and directly imaged gas giant exoplanet atmospheres and clouds. Atmospheric structure can be investigated by monitoring ultracool object rotational modulation (Buenzli et al. 2012, 2015a, 2015b; Apai et al. 2013; Karalidi et al. 2016; Yang et al. 2016; Schlawin et al. 2017; Biller et al. 2018; Zhou et al. 2020a; Manjavacas et al. 2021). The Hubble Space Telescope (HST) Cloud Atlas program has focused on using these methods and performed extensive spectrophotometric studies of BDs, gas giants, and planetary-mass objects (Lew et al. 2016, 2020b, 2020a; Manjavacas et al. 2018, 2019b, 2019a; Zhou et al. 2018, 2019, 2020b; Miles-Páez et al. 2019). Moving forward, these methods can potentially be integrated with techniques such as Doppler imaging to create more complete surface maps for ultracool objects. The unified model presented in this paper seeks to address this need by combining photometric techniques with Doppler imaging data sets.

In this paper, we develop analytical and numerical frameworks for Doppler imaging ultracool objects using the Euler−Rodrigues formula (Section 2). We internally validate our analytical method versus a numerical stellar/substellar surface model (Section 3). Our model is then externally validated against spectroscopic and photometric archival data from Crossfield et al. (2014) and Buenzli et al. (2015b) (Section 4). Based on our results, we discuss the implications for Doppler imaging of ultracool objects and future avenues of research (Section 5).

LPs have a higher S/N than individual spectral lines and are therefore more effective in identifying surface inhomogeneities. One common method of computing LPs for spectra uses the cross-correlation function; however, for fast-rotating objects, least squares deconvolution (LSD) can be a mathematically rigorous technique. LSD is an attractive option for both spectroscopic and spectropolarimetric applications because it can create high-precision, averaged LPs, ideal for Bayesian inference and Doppler imaging (Kochukhov et al. 2010).

To understand how LSD works, we must consider how the target's observed spectral lines are broadened through multiple processes including rotation and instrumentation. We can account for these processes mathematically by convolving the source spectrum S with a broadening and instrument kernel, represented by Z. This results in the observed spectrum Y⁰ determined as such,

$$\begin{eqnarray}&&{Y}^{0}=S\ast Z,\end{eqnarray} \tag{ 1 }$$

where ∗ is the convolution operator. In Equation (1), the convolution operator broadens the source spectrum, S, by the width of the broadening kernel, Z, creating a mathematical representation of the observed spectrum.

In most cases, we do not have perfect information and cannot directly observe the broadening kernel. The source spectrum is only estimated based on prior observations and models. The parameters we wish to determine are often innate to the source spectrum, but our primary observable is the observed spectrum, Y⁰. We estimate the source spectrum with a template spectrum, F, which can be computed from various models based on known or estimated parameters such as effective temperature (T_eff), log(g), and metallicity.

We conducted sensitivity tests for both effective temperature and log(g) with a fixed template spectra (T_eff = 1500 K, log(g) = 5.0, [Fe/ H] = 0), in the range used in Crossfield et al. (2014). We created synthetic observed spectra with varying effective temperature (±50 K) and log(g) (±0.5); available spectra in this temperature regime did not have multiple metallicities available. This synthetic spectra was implemented, along with the fixed template, into our numerical model (see Section 2.3) to create LPs with a spot of known latitude, longitude, contrast, and radius. We then retrieved these parameters via Bayesian inference. A temperature difference of 50 K resulted in the following deviations in the retrieved parameters: latitude (1.24σ), longitude (0.57σ), radius (0.57σ), and contrast (0.19σ). A log(g) difference of 0.5 resulted in the following: latitude (1.68σ), longitude (0.096σ), radius (0.9σ), and contrast (0.16σ). As can be seen, longitude, radius, and especially contrast were fairly robust to temperature and gravity deviations while latitude was less robust.

We can replace S with F in Equation (1) as follows:

$$\begin{eqnarray}&&{Y}^{0}=F\ast Z={MZ}.\end{eqnarray} \tag{ 2 }$$

In this equation, Y⁰ is comprised of m elements with spacing dependent on the spectral sampling of the observational instrument. In practice, solving Equation (2) requires expanding F into a line mask, M, via an m × n Toeplitz matrix as described in Donati et al. (1997). This provides us with the following form:

$$\begin{eqnarray}&&Z={M}^{-1}\ {Y}^{0}.\end{eqnarray} \tag{ 3 }$$

Similar to Crossfield et al. (2014), we use LSD to convert the observed spectrum into an LP in RV space. These LPs can then be compared to either a mean LP from the observation or to a modeled LP to infer surface inhomogeneities. The specific LSD method used in this paper was adopted from Wang et al. (2017, 2018) and Pai Asnodkar et al. (2022), and is best represented by the following:

$$\begin{eqnarray}&&Z={\left({M}^{T}\ {S}^{2}\ M+{\rm{\Lambda }}R\right)}^{-1}\ {M}^{T}\ {S}^{2}\ {Y}^{0},\end{eqnarray} \tag{ 4 }$$

where Z is the computed LP comprised of n data points and S is an m × m diagonal error matrix with S_ii elements equal to 1/σ_ii where σ_ii is the measurement error at each epoch. Λ is a scalar regularization parameter while R is the associated regularization Toeplitz matrix, which is outlined in Donatelli & Reichel (2014).

2.2.1. Scheme

Doppler tomography examines the changes in a target's LP during a transit of its viewable hemisphere. Our analytical method for modeling stellar and substellar surfaces derives inspiration from Doppler tomographic studies of eclipsing binary systems (Albrecht et al. 2007; Wang et al. 2018) and transiting exoplanets (Hirano et al. 2011; Johnson et al. 2014). Spectral deviations created by a transiting companion or surface inhomogeneity blocking a portion of the target's visible surface can be measured using a high-S/N technique such as LSD.

These spectral deviations create Doppler shadows when they are plotted in RV space. With this technique, companion stars and planets will create linear shadows when they transit the target. In this article, we make the unique contribution of considering that surface inhomogeneities create curved Doppler shadows, if the target is inclined with respect to the viewer, as can be seen in Figure 1. For the hemisphere inclined toward the viewer, higher latitudes spots are more extensive in their viewable phase.

Zoom In Zoom Out Reset image size

We represent the unperturbed stellar/substellar disk as a rotationally broadened LP. Surface inhomogeneities are modeled as Gaussian distributions whose width depends on radius and viewing angle. These Gaussian distributions are either subtracted from (dark spot) or added to (bright spot) the rotationally broadened LP to create an analytical model for each epoch. An example of this process can be seen in Figure 2 with dark spots subtracting from the rotationally broadened LP.

Zoom In Zoom Out Reset image size

2.2.2. Unperturbed LP

We approximate the unperturbed LP with a rotational broadening kernel. Rotation is the dominant spectral line broadening mechanism for fast-rotating astrophysical objects ($\upsilon \ \sin \ i\gt 20\ \mathrm{km}\,{{\rm{s}}}^{-1}$, where υ is the equatorial rotation speed of the object and i is its inclination).

A suitable expression for this rotational broadening kernel is described in Gray (2008) and adapted by Pai Asnodkar et al. (2022) as follows:

$$\begin{eqnarray}&&\begin{array}{l}G({\upsilon }_{i},t)\\ ={c}_{1}{\left[1-{\left(\displaystyle \frac{{\upsilon }_{i}-{\upsilon }_{\star }(t)}{\upsilon \ \sin \ i}\right)}^{2}\right]}^{1/2}\\ +\,{c}_{2}\left[1-{\left(\displaystyle \frac{{\upsilon }_{i}-{\upsilon }_{\star }(t)}{\upsilon \ \sin \ i}\right)}^{2}\right].\end{array}\end{eqnarray} \tag{ 5 }$$

G(υ_i , t) is the rotational broadening kernel applied over velocities υ_i and epochs, t, while υ_⋆ is the RV of the star with respect to the solar system. The constants c₁ and c₂ are computed as follows:

$$\begin{eqnarray}&&{c}_{1}=\displaystyle \frac{2\left(1-\epsilon \right)}{\pi \upsilon \ \sin \ i\left(1-\epsilon /3\right)}\end{eqnarray} \tag{ 6 }$$

and

$$\begin{eqnarray}&&{c}_{2}=\displaystyle \frac{\epsilon }{2\upsilon \ \sin \ i\left(1-\epsilon /3\right)},\end{eqnarray} \tag{ 7 }$$

where is the linear limb darkening coefficient describing how the center of the stellar or substellar target appears brighter than the edge.

As the linear limb darkening law represents a poor approximation of the limb darkening phenomenon for stellar and substellar targets, we conducted a comparison between the rotational broadening kernel in Equation (5) and one created via higher-order Claret limb darkening. Using Scipy's least squares tool (Virtanen et al. 2020), we fit Equation (5) to an LP created with Claret's limb darkening coefficients for a 1500 K, log(g) = 5.0, solar metallicity, substellar object (Claret et al. 2012). Both profiles were also convolved with an instrument profile. The resulting profiles agreed within 0.2%, demonstrating the suitability of Equation (5) for representing ultracool stellar and substellar disks.

2.2.3. Gaussian Spots

Surface inhomogeneities can be described by Gaussian distributions centered on an RV value that is a function of the longitudinal viewing angle with respect to the observer. The Gaussian distribution's width is set by the surface inhomogeneity's longitudinal radius. The magnitude of the Gaussian distribution is a product of its area and contrast.

As the target star or substellar object rotates, the Gaussian distribution will transit across the LP. The range of the transit in RV space depends on the inhomogeneity's latitude. As demonstrated in Figure 1, an equatorial inhomogeneity will transit the entire rotationally broadened LP while an inhomogeneity at mid-to-polar latitudes will transit an RV band commensurate with its latitude. Conversely, a polar spot on an inclined target may not rotate out of view such that while its RV domain is limited by its latitude, its signature will be viewable throughout the entire phase space.

We follow the examples of Crossfield et al. (2014), Karalidi et al. (2015), and Karalidi et al. (2016), by modeling surface inhomogeneities as circles or ellipses. These inhomogeneities are initially modeled as two-dimensional Gaussian ellipses with the longitudinal and latitudinal widths adjusted for viewing angle by Lambert's cosine law,

$$\begin{eqnarray}&&I(\theta )={I}_{0}\ \cos (\theta ),\end{eqnarray} \tag{ 8 }$$

which states that the observed intensity of an Lambertian object decays according to a cosine function where θ is the angle of the observer with respect to the normal. Assuming a circular Gaussian spot, the longitudinal (r_L ) and latitudinal (r_l ) radii then become

$$\begin{eqnarray}\begin{array}{rcl}{r}_{L} & = & r\ \cos (\alpha )\\ {r}_{l} & = & r\ \cos (\beta ),\end{array}\end{eqnarray} \tag{ 9 }$$

where r is the radius of the Gaussian spot with units of degrees and α and β are the longitudinal and latitudinal viewing angles. These viewing angles are measured from the sub-observer centroid of target to the center of the Gaussian spot. They are computed by first assuming that the target is a sphere with radius equal to unity and the following spherical coordinates:

$$\begin{eqnarray}\begin{array}{rcl}r & = & 1\\ \theta & = & 90-l\\ \phi & = & L.\end{array}\end{eqnarray} \tag{ 10 }$$

The spot's latitude and longitude are then converted into Euler coordinates, $\vec{x}=(x,y,z)$,

$$\begin{eqnarray}\begin{array}{rcl}x & = & r\ \cos (\phi )\ \sin (\theta )\\ y & = & r\ \sin (\phi )\ \sin (\theta )\\ z & = & r\ \cos (\theta ).\end{array}\end{eqnarray} \tag{ 11 }$$

We define the rotational angle (ξ) as a function of the inclination (ξ = 90°–i). The transformation matrix (T) is composed of rotation parameters (a, b, c, d). Computing these rotation parameters, we encounter the components of the unit vector $\vec{k}=({k}_{x},{k}_{y},{k}_{z})$, but as we will only rotate along the y-axis by varying the inclination, k_x = k_z = 0 (Shuster 1993).

$$\begin{eqnarray}\begin{array}{rcl}a & = & \cos \left(\displaystyle \frac{\xi }{2}\right)\\ b & = & {k}_{x}\ \sin \left(\displaystyle \frac{\xi }{2}\right)=0\\ c & = & {k}_{y}\ \sin \left(\displaystyle \frac{\xi }{2}\right)\\ d & = & {k}_{z}\ \sin \left(\displaystyle \frac{\xi }{2}\right)=0.\end{array}\end{eqnarray} \tag{ 12 }$$

The rotated Euler coordinates $(\vec{x^{\prime} })$ can now be computed with the following transformation (Shuster 1993):

$$\begin{eqnarray}&&\vec{x^{\prime} }=T\vec{x},\end{eqnarray} \tag{ 13 }$$

where T is defined as follows (Shuster 1993):

$$\begin{eqnarray}&&\begin{array}{l}T=\\ \left[\begin{array}{ccc}{a}^{2}+{b}^{2}-{c}^{2}-{d}^{2} & 2({bc}-{ad}) & 2({bd}+{ac})\\ 2({bc}+{ad}) & {a}^{2}+{c}^{2}-{b}^{2}-{d}^{2} & 2({cd}-{ab})\\ 2({bd}-{ac}) & 2({cd}+{ab}) & {a}^{2}+{d}^{2}-{b}^{2}-{c}^{2}\end{array}\right].\end{array}\end{eqnarray} \tag{ 14 }$$

The Euler−Rodrigues formula is used to account for the coordinate transformation necessitated by any orbital inclination that is not 90°. The transformation results in new coordinates ($\vec{x^{\prime} }$), which can then be converted back into the longitudinal and latitudinal viewing angles, α and β. Using Equation (9), the longitudinal radius, r_L , is then converted into RV space for use in the LP; this will be the standard deviation (σ) of the Gaussian spot.

Using the Euler−Rodrigues formula to compute a spot's position with respect to the observer's viewing angle creates a set of coordinates which, depending on the duration of observations, could include spot positions encompassing an entire revolution. In subsequent analysis, we will compute the longitudinal and latitudinal viewing angles for one complete revolution. Depending on the target's inclination, along with the spot's latitude, the spot may have instances where it is not viewable for periods of time. In the opposite case, for positive inclinations (0° < i < 90°), if the latitude is greater than inclination (l > i), spots in the Northern Hemisphere (hemisphere tilted toward the observer) will remain in view for the entire revolution as they circle the North Pole. Each of these cases must be taken into account and spot-induced LP deviations will be set to zero when the spot is not within view of the observer.

The magnitude of the Gaussian spot is a function of the spot's area and contrast. The area (a_spot) is computed as a fraction of the star's unit disk:

$$\begin{eqnarray}&&{a}_{\mathrm{spot}}=\displaystyle \frac{{r}_{L}\ {r}_{l}}{{\left(90^\circ \right)}^{2}}.\end{eqnarray} \tag{ 15 }$$

This leads us to the Gaussian distribution for the surface inhomogeneity:

$$\begin{eqnarray}&&p({RV})=\displaystyle \frac{{a}_{{\rm{spot}}}\cdot {\rm{contrast}}}{\sigma \sqrt{2\pi }}\ {e}^{-\tfrac{1}{2}{\left(\tfrac{{RV}-{{RV}}_{{\rm{spot}}}}{\sigma }\right)}^{2}},\end{eqnarray} \tag{ 16 }$$

where

$$\begin{eqnarray}&&\sigma ={r}_{L}\displaystyle \frac{v\ \sin \ i}{90^\circ }\end{eqnarray} \tag{ 17 }$$

and

$$\begin{eqnarray}&&{{RV}}_{\mathrm{spot}}=\alpha \displaystyle \frac{v\ \sin \ i}{90^\circ }.\end{eqnarray} \tag{ 18 }$$

RV_spot is the spot centroid RV for each epoch and RV ranges from −3σ to 3σ to capture the continuum well beyond the expected rotationally broadened signal at ±1σ for normalization purposes. The term v sin i/90° is used to convert from longitude to velocity space.

Contrast varies in the range [−1, 1] with −1 being twice as bright as the average surface, 0 demarcating the average surface temperature, and 1 denoting a perfectly dark spot with 0 K temperature and no photon emission. It can be seen that a negative contrast will lead to more flux added to an LP centered at the spot's location while positive contrast will lead to a subtraction of flux centered at the spot's location. The spot's maximum achievable magnitude is simply the product of its area and its contrast.

For each epoch, the computed Gaussian spot in the RV domain is then added or subtracted to the rotational broadening LP. Figure 2 demonstrates the analytical method with a starspot moving across a 60° inclined disk with an equatorial latitude; −45°,0°, and +45° longitudes; 10° radius; and 0.5 contrast.

2.2.4. Analytical Model Photometry

While an analytical starspot model has clear advantages in terms of elegance and computational speed, a numerical model for stellar and substellar photospheric surfaces is necessary to verify the validity of the analytical approach. Building a numerical model for the star also provides a convenient option for creating global maps of the star's surface to verify each spot's location as well as their motion (ensuring that the motion computed by the analytical method is indeed physical and passes the so-called visual check).

To create the numerical model, the star's surface is divided n times longitudinally and n/2 times latitudinally to create n²/2 cells. To speed computation, n longitudinal cells are implemented so that each longitudinal cell directly maps to an RV value. For example, a longitude of +90° directly maps to an RV of + υsin i. Likewise, a longitude of −90° maps to an RV of −υ sin i.

It is also necessary to consider the unseen hemisphere. The observer only has a view of one hemisphere, so the flux at any point located at an angular distance greater than 90° from the sub-observer point (α = 0°, β = 0°), must be set to zero.

All the points within the 90° radius are initially set to 1 before being multiplied by a linear limb darkening law (to align with the linear limb darkening used in the rotational broadening equation in Section 2.2 as well as to match the practice in Crossfield et al. 2014). The insets of Figures 1 and 2 demonstrate the numerical code's ability to produce limb-darkened stellar surfaces.

Star spots can then be added to the numerical model using a bivariate Gaussian distribution with the longitudinal and latitudinal radii (standard deviations) calculated with Equation (9). As in the analytical technique, the Euler−Rodrigues formula is used to transform the spot's coordinates based on the target's inclination for one entire rotation. The spots are then added or subtracted (bright spots are added while dark spots are subtracted) to the limb-darkened image.

The numerical code can be implemented to generate a dynamic LP that accounts for the rotation of spots into and out of view as well as the target's inclination. As mentioned above, the target's longitude directly maps to the RV in the LP. Figure 3 demonstrates the starspot's effect on the LP. Both the numerical and analytical codes lead to a smaller absorption profile in the redshifted hemisphere, which aligns with our expectations that the LP should experience a blueshift when a surface inhomogeneity is in the Eastern Hemisphere (the hemisphere rotating away from the observer).

Zoom In Zoom Out Reset image size

Similar to the analytical model, the numerical code can also be used to generate photometric light curves. This is accomplished by summing up the total flux across the numerical map at each epoch. This will create a value that far exceeds unity, so normalization is necessary. To accomplish this normalization, a featureless numerical map of the star or substellar object is generated using the selected limb darkening law. The sum of flux at each epoch is then divided by this featureless sum of flux to create a normalized value for each measurement. A light curve closely matching the analytical photometry is produced, as seen in Figure 4.

Zoom In Zoom Out Reset image size

A similar issue arises with the numerical curve as that seen with the analytical curve, the spot's curve is not smooth. We apply the same solution and use a 1D Gaussian filter to smooth the curve leading to a more physical appearing light curve that also aids in fitting data with Bayesian inference.

We use the open-source Python module, Dynesty, to infer spots' locations, size, and brightness using our analytical model as well as assumed priors (Speagle 2020). For our purposes, performing Bayesian inference on stellar/substellar surfaces using Markov Chain Monte Carlo (MCMC) methods does not necessarily result in optimal or easy-to-decipher posterior probabilities due to the multimodal nature of star spots and BD atmospheric features once one moves past a 1-Spot model. For this reason, a nested sampling approach was adopted. Dynesty not only functions well with multimodal data, it also computes the marginal likelihood, or in nested sampling terms, the evidence (Speagle 2020).

John Skilling developed nested sampling to solve for the evidence by computing the likelihood and prior mass volume through series of nested volumes of likelihood (Skilling 2004, 2006). Live points are distributed throughout the prior volume and at each iteration the live point with the lowest likelihood is replaced with a new live point with a greater likelihood; the discarded live points (dead points) are then used to compute the evidence and posterior distribution (Skilling 2004, 2006).

In this section, we will perform a sequence of validation tests, each resolving two spots with varying angular separation and inclination.

3.2.1. Test Cases Based on Karalidi et al. (2015)

To accomplish a systematic exploration of the model's robustness, we adopt a similar approach to that of Karalidi et al. (2015), which validated Aeolus, a photometric MCMC-based substellar surface mapping code. Karalidi et al. (2015) tested their code by varying the angular distance in terms of latitude and longitude for two spots with an approximately infinite S/N. Spot coordinates were input and the model was used to retrieve both spots location. For this work, the numerical model was implemented with input spot locations, sizes, and contrasts to generate simulated observed LPs. We then performed a Bayesian inference using Dynesty to determine the spot locations.

As in Karalidi et al. (2015), we assume a larger spot (fixed at the equator at the observer's latitude) with a radius of 9° and contrast of 0.3. The smaller spot (with variable location) has a radius of 5° and contrast of 0.6. For this validation test, $v\ \sin \ i=50\ \mathrm{km}\,{{\rm{s}}}^{-1}$ was selected based on the rotational range recommended by Vogt et al. (1987), and the inclination is nearly edge-on with i = 85°.

For each case, the input locations are retrieved within 1σ for both spots. As can be seen in Figure 5, the analytical model in this work improves on the angular distance that can be resolved: both the smaller angular separations of 24° longitudinally and 20° latitudinally can be resolved.

Zoom In Zoom Out Reset image size

We note that the reader may notice in Figure 5 the angular distances from Karalidi et al. (2015) do not appear to exactly match the distances from the analytical model from this work. This is an artifact relating to how the maps in Karalidi et al. (2015) are constructed with two-dimensional circles while the maps generated in this work are effectively three dimensional due to the Euler−Rodrigues transformation. For example, an equatorial spot traveling 10° across a stellar disk will appear to travel less because of the larger radial motion along the line of sight than the tangential motion.

Another important difference between the two models is this paper's adoption of Lambert's cosine law. In Figure 5, the effects of the law can be seen in how the spot's shape is modified near the stellar limb to account for the foreshortening of the received flux for a spot seen at high viewing angles.

An example of the input simulated LP data and the retrieved 2-Spot analytical model for the ΔL = 24° case can be seen in Figure 6. Deviations for both spots can be seen to transit the LP during the early and late phases of the rotation. The 1% noise level can also be observed in the simulated data.

Zoom In Zoom Out Reset image size

3.2.2. Comprehensive Grid Test

To truly replicate results in Karalidi et al. (2015), we needed to compare our photometric model to that shown in their work. To this end, we performed the same analysis as shown in Section 3.2 only with photometric data.

Again, we set the rotational speed to v sin i to 50 km s⁻¹ and the viewing angle to nearly edge-on, i = 85°. The primary and secondary spots were at the same locations with identical spot radii and contrasts as outlined in Section 3.2. The numerical model was used to generate a photometric light curve as described in Section 2.3 while the analytical model was used to compute the likelihood within a nested sampling Bayesian inference. Both the analytical and numerical photometric light curves were smoothed using a Gaussian filter as described in Sections 2.2 and 2.3.

The results are very similar to those shown in Figure 5 with improved spot resolution versus that seen in Karalidi et al. (2015). Within the context of this sample case, the analytical code resolves both spots' locations for both spectroscopic and photometric data.

As in Section 3.2, we conducted an exhaustive exploration for varying longitude and latitude displacements for i = 30° and 60°. The spot size and contrast as well as the noise level (1%) were identical to those discussed in Section 3.2.2. Just as in that section, the photometric analytical model successfully retrieved every input scenario.

Luhman 16AB are two of the solar system's nearest neighbors at 1.9955 ± 0.0004 pc, which also makes the two binary components the nearest BDs (Bedin et al. 2017). Luhman 16AB was originally discovered by Luhman (2013) using images from the Digitized Sky Survey, the Two Micron All-Sky Survey, and the Deep Near-Infrared Survey of the Southern Sky. Luhman 16AB is comprised of Luhman 16A, a ${34.2}_{-1.1}^{+1.3}\ {M}_{J}$ primary with an L7.5 spectral class, and Luhman 16B, a ${27.9}_{-1.1}^{+1.0}{M}_{J}$ secondary with a measured spectral class T0.5 (Burgasser et al. 2013; Kniazev et al. 2013; Garcia et al. 2017; Ammons & Garcia 2019).

Similar to many BDs at the L/T transition, Luhman 16B demonstrates high photometric and spectroscopic variability (Biller et al. 2013; Gillon et al. 2013; Burgasser et al. 2014; Buenzli et al. 2015a, 2015b). Although the primary, Luhman 16A, is also considered variable, it is to a smaller degree than the secondary, Luhman 16B (Biller et al. 2013; Burgasser et al. 2014; Crossfield et al. 2014; Buenzli et al. 2015a). Several studies over the last decade have investigated Luhman 16B's variability and periodicity. Gillon et al. (2013) deduced an 11% peak-to-peak near-infrared (0.9 μm) variability and 4.87 ± 0.01 hr period over a 12 night campaign. Burgasser et al. (2014) estimated a similar 13.5% near-infrared (1.25 μm) variability. Apai et al. (2021) used TESS to observe 100 rotations of the binary system, finding maximum peak-to-peak variability to be approximately 4% with longer period, 10% variability seen over 100 hr timelines. Apai et al. (2021) attributed the strongest periodicity signal (5.28 hr) along with a secondary periodicity (2.5 hr) to the combined effects of differential rotation and planetary waves in Luhman 16B's atmosphere. In line with the above results, Heinze et al. (2021) observed a 5 hr period and variability in the red continuum changing from 7% to 13% on consecutive nights, providing a relatively consistent measured period and high-variability.

Global surface maps of Luhman 16B's surface inhomogeneities have been generated from both spectroscopic and photometric observations. Crossfield et al. (2014) inferred spot latitude, contrast, and size based on a 1-Spot model. They also created the first surface map of the BD by applying the maximum entropy method outlined in Vogt et al. (1987).

More recently, Luger et al. (2021) reprocessed the Crossfield et al. (2014) data and used the open-source Python code Starry (Luger et al. 2019) to generate a global surface map. The map created by Luger et al. (2021) is qualitatively similar to that of Crossfield et al. (2014) with both studies deducing an dark equatorial region and a bright polar spot.

Using photometric data from the HST, Karalidi et al. (2016) generated three and four spot surface maps for Luhman 16A and B. For Luhman 16B, a four spot model with temperature (contrast) as a free parameter was qualitatively different from the maps created by Crossfield et al. (2014) and Luger et al. (2021), but still contained near equatorial spots and brighter spots in the polar region (Karalidi et al. 2016). In the following subsections, we will use archival data from both of these studies to produce Luhman 16B surface maps with the analytical method described in this paper.

To infer a global map of Luhman 16B's surface, we follow a process similar to that used by Crossfield et al. (2014) with the same VLT/CRIRES data collected on 2013 May 5. Crossfield et al. (2014) initially obtained 56 spectra (2.288−2.345 μm) with 300 s exposure times and combined these into groups of four to create high S/N LPs via the LSD method. We obtained the Luhman 16B data via open online archival and used the reprocessed template and observed spectra from Luger et al. (2021) to similarly create LPs via the process outlined in Section 2.1. Referencing Section 2.1, the template spectrum created by Crossfield et al. (2014) and reprocessed by Luger et al. (2021) is used to create a line mask matrix (M). Equation (4) is then used to compute the LP for each epoch.

An important point to note is the retrieved latitude's sensitivity to factors such as LP element number (n) and the regularization parameter (Λ). We selected n = 93 LP elements to ensure a ΔRV ≈ 1−2 km s⁻¹ while also sampling the RV range of ±3σ. Λ = 500 was qualitatively selected to attenuate noise in the retrieved LP while preserving the observed signal. Varying these two parameters can vary the retrieved latitude, but all resulting variations of the retrieved latitudes are consistent with those of Crossfield et al. (2014).

The resulting LPs can be seen in Figure 7 along with the line fit from our 1-Spot analytical model. Due to the retrieved spot's radius and relatively low contrast (compared to the validation examples), the affected region of the LP is broad and the subtle shifts in flux are difficult to detect by eye, necessitating a computational Bayesian method. Using nested sampling, we inferred spot longitude, latitude, radius, and contrast by comparing the LPs generated from observed spectra to our analytical model. As in Crossfield et al. (2014) and Luger et al. (2021), we set inclination to i = 70° and v sin i = 26.1 km s⁻¹. We manually varied the linear limb darkening coefficient () and the number of LP elements (n) to fit the observed data. We convolved the LP with a Gaussian profile to account for instrument and other sources of line broadening resulting in the modeled LPs in Figure 7.

Zoom In Zoom Out Reset image size

For the 1-Spot model, dynamic nested sampling resulted in an equatorial dark spot (latitude = 29°) with radius = 37° and contrast = 0.10, which can be seen in Figure 8. These results can be compared to those of Crossfield et al. (2014) and Karalidi et al. (2016), and the photometric results in Table 1. Our results match well with the values given in Crossfield et al. (2014), which inferred a latitude ≤31°, radius = 33° ± 7°, and contrast = 0.12 ± 0.03. Longitude was not included in Crossfield et al. (2014) as it is dependent on the beginning epoch of observation.

Zoom In Zoom Out Reset image size

Table 1. Summary of the Model Comparison ^a

Model	Mode	Latitude	Longitude	Radius	Contrast
Crossfield et al. (2014)	Spectroscopic	≤31°	N/A	33° ± 7°	0.12 ± 0.03
1-Spot	Spectroscopic	$29^\circ \frac{+8.44^\circ }{-8.90^\circ }$	$-160.39^\circ \frac{+2.08^\circ }{-2.21^\circ }$	$36.79{^\circ }_{-6.73^\circ }^{+7.73^\circ }$	${0.10}_{-0.03}^{+0.04}$
Karalidi et al. (2016)	Photometric (J)	−20° ± 12°	113.47° ± 6.46°	38.89° ± 0.8°	−0.180 ± 0.013
	(fixed bright spots)	45° ± 12°	186.08° ± 3.65°	24.45° ± 1.46°	−0.180 ± 0.013
		72° ± 10°	283° ± 13°	35.54° ± 1.67°	−0.180 ± 0.013
Karalidi et al. (2016)	Photometric (H)	−18° ± 8°	107.27° ± 5.48°	32.42° ± 0.8°	−0.180 ± 0.013
	(fixed bright spots)	51° ± 20°	190.91° ± 4.83°	25.02° ± 1.13°	−0.180 ± 0.013
		74° ± 6°	288.6° ± 9.4°	36.30° ± 1.18°	−0.180 ± 0.013
1-Spot	Photometric (J/H)	$-23.53^\circ \frac{+15.72^\circ }{-9.89^\circ }$	$110.08{^\circ }_{-2.79^\circ }^{+3.30^\circ }$	$37.58{^\circ }_{-5.37^\circ }^{+4.86^\circ }$	${0.25}_{-0.06}^{+0.07}$
2-Spot	Photometric (J/H)	$-27.20{^\circ }_{-10.86^\circ }^{+18.22^\circ }$	$119.30{^\circ }_{-5.22^\circ }^{+5.92^\circ }$	$34.86{^\circ }_{-4.85^\circ }^{+3.59^\circ }$	${0.28}_{-0.06}^{+0.05}$
		$47.19{^\circ }_{-32.52^\circ }^{+24.41^\circ }$	$\ast 231.68{^\circ }_{-20.44^\circ }^{+15.54^\circ }$	$28.20{^\circ }_{-7.44^\circ }^{+7.64^\circ }$	$-{0.14}_{-0.09}^{+0.06}$

^a Note. ∗ Denotes where +360° has been added to negative longitudes to more easily compare to external results where the longitudinal phase goes from 0° to 360° vs. −180° to 180°

Download table as:  ASCIITypeset image

The 1-Spot model reproduces the MCMC results in Crossfield et al. (2014) and shows them to be a reasonable representation of Luhman 16B's dominant feature, the equatorial dark spot (Figure 9). For models with a higher number of spots, the spectroscopic analytical method did not converge on a likelihood, leading us to believe the 1-Spot model best describes the Crossfield et al. (2014) data.

Zoom In Zoom Out Reset image size

Luhman 16A and B photometric data were collected on 2013 November 8 by the HST/Wide-Field Camera with the G141 grism, which operates between 1.1 and 1.66 μm. The data was originally collected, processed, and reported by Buenzli et al. (2015b). Karalidi et al. (2016) used those observations, as well as those collected by Buenzli et al. (2015a) with the G102 grism (0.8−1.15 μm) on 2014 November 23, to construct three and four spot surface maps using the Aeolus code developed in Karalidi et al. (2015).

Here we will seek to infer surface spots from the Buenzli et al. (2015b) data using the analytical method described in Section 2.2. The J and H near-infrared bands from Buenzli et al. (2015b) result in very similar light curves, so for this paper these bands are averaged.

In Figure 10, the Buenzli et al. (2015b) J- and H-band data, along with our 1- and 2-Spot photometric models are plotted. The 2-Spot light-curve model more fully represents Luhman 16B's photometric behavior. Dynamic nested sampling results in a model with improved evidence (${\rm{\Delta }}\mathrm{ln}{ \mathcal Z }=13.8$) for the 2-Spot model over the 1-Spot model. This is strong evidence that the 2-Spot model better represents the data (Kass & Raftery 1995). A summary of the photometric results can be found in Table 1.

Zoom In Zoom Out Reset image size

Further increasing the number of spots in our model leads to lower comparative evidence (${\rm{\Delta }}\mathrm{ln}{ \mathcal Z }$), so we conclude the strongest evidence exists for the 2-Spot model. The equatorial dark spot and midlatitude bright spot inferred by the 2-Spot photometric model echo the results seen in the photometric results in Karalidi et al. (2016) as well as the Doppler imaging maps in Crossfield et al. (2014) and Luger et al. (2021).

A summary of both spectroscopic and photometric results, as well as comparisons to those of Crossfield et al. (2014) and Karalidi et al. (2016), can be seen in Table 1. Both the spectroscopic and photometric results reasonably match those obtained in Crossfield et al. (2014) and Karalidi et al. (2016), respectively. Even though the results from Karalidi et al. (2016) rely on fixed bright spots (negative contrast for this paper's nomenclature), our 1-Spot and 2-Spot analytical models still approximately infer the average spot location and size of the J and H bands in Karalidi et al. (2016). This matches our expectations as we averaged the J and H bands before performing Bayesian inference via nested sampling. Furthermore, our 2-Spot photometric model resulted in a Bayesian information criterion value of 41.4, matching the J-band 2-Spot model value achieved in Karalidi et al. (2016).

Multiple data sets at different wavelengths have explored the surface of Luhman 16B. From these studies, evidence has emerged that the BD possesses a persistent dark equatorial spot and potentially a bright mid-to-polar latitude spot.

Persistent atmospheric features for Luhman 16B can be seen in the results of this paper as well as those in Crossfield et al. (2014), Luger et al. (2021), and Karalidi et al. (2016). As discussed above, Crossfield et al. (2014) inferred a dark equatorial spot using MCMC. Using the maximum entropy method, Crossfield et al. (2014) also created a global map that included a dark equatorial spot as well as a bright polar feature. Luger et al. (2021) reprocessed the same spectroscopic data and inferred a similar map with those same features. Karalidi et al. (2016) identified the dark feature and labeled it the Possible Persistent Cloud Structure (PPCS-1), which may be similar to Jupiter's Great Red Spot. Not only did it appear in the HST J and H bands as well as a G102 grism light curve, it also appears in data presented in Gillon et al. (2013), which used the I + z TRAPPIST filter (Karalidi et al. 2016).

These features may be consistent with those in Apai et al. (2021), in which Luhman 16B's planetary waves, potentially in the form of zonal jets, were demonstrated over 100 hr of TESS observations. It is possible the persistent dark spot forms in such a zonal jet at Luhman 16B's equator while mid-to-polar spots form within higher latitude bands. Apai et al. (2021) also detected what they thought to be a long period signal from a notional polar vortex on Luhman 16A. The existence of such persistent features could help shed light on the overall atmospheric structure and dynamics of Luhman 16B and similar ultracool objects.

While there is evidence for persistent features, variability remains a recurring theme in observations of BDs. L/T transition variability is likely caused by multiple atmospheric layers and clouds with disparate opacities. The variability of Luhman 16B (Gillon et al. 2013; Crossfield et al. 2014; Buenzli et al. 2015a, 2015b) may also arise from the different observational techniques and wavelengths considered. Not only did these studies observe Luhman 16B at different times, but they also probed different pressure levels and chemical species. The rich inhomogeneous data sets can complicate the picture when attempting to quantify the true variability of the object. However, they also offer an opportunity to understand the BD's three-dimensional atmospheric dynamics.

Over the past decade, there has been extensive work in modeling atmospheric dynamics in lightly irradiated BD atmospheres through general circulation models (GCMs). Several drivers for robust circulation in the stratified overlying atmosphere of these BDs have been explored in the literature, including convection driving the stratified upper layer and cloud-radiative feedback. Showman et al. (2020) and Zhang (2020) both include comprehensive reviews of these efforts.

Convection from the interior perturbing the overlying stratified atmospheric layer has been a mechanism identified as a potential primary driver of BD circulation (Freytag et al. 2010; Showman & Kaspi 2013). The strength of these perturbations, along with the object's radiative timescale and drag at the radiative-convective boundary, appear to control the latitudinal extent and strength of jet stream formation (Zhang 2014; Showman et al. 2019; Tan 2022).

Cloud-radiative feedback has also been demonstrated as a plausible mechanism for driving atmospheric circulation. One-dimensional feedback between clouds, radiation, and atmospheric mixing is explored in Tan (2019). More recently, modelers have focused on three-dimensional feedback models in which cloud-radiative dynamics induce active atmospheric circulation thereby perpetuating cloud patchiness (Tan 2021a, 2021b; Lefèvre et al. 2022).

Placing Luhman 16B's large equatorial dark spot and midlatitude/polar bright spots within the context of a model such as that discussed in Zhang (2014) offers a means to connect Doppler imaging to atmospheric dynamics. In that work, the atmospheric structure is essentially a product of forcing from the convective interior and the radiative timescale of the overlying stratified atmosphere. Strong convective forcing and weak radiation lead to large equatorial jets with midlatitude to polar vortices and turbulence. It is possible that the equatorial dark spots are persistent clouds or storm systems within a jet stream band. The bright spots may be regions of lessened opacity and cloud thickness created by turbulent features in which we can observe deeper layers of the BD's atmosphere.

Combining BD atmospheric circulation models with Doppler imaging techniques could potentially lead to more accurate surface maps as well as validation for the development and selection of GCMs. Follow-on observations of Luhman 16B and other ultracool targets may result in more accurate atmospheric models.

Doppler imaging focused on specific chemical species' spectral lines could help to further illuminate the mechanisms at play in BD and gas giant exoplanet atmospheres. Tracing specific species such as CO, CH₄, NH₃, and molecular nitrogen (N₂), provides insight into atmospheric circulation. At temperatures above the L/T transition, carbon's chemical equilibrium state transitions from favoring CO (above 1470 K) to CH₄ (below 1470 K) (Fegley & Lodders 1996). A similar chemical equilibrium transition occurs for nitrogen. Below 630 K, NH₃ should be the dominant species while N₂ is favored by chemical equilibrium above this value (Fegley & Lodders 1996).

Chemical disequilibrium in ultracool atmospheres provides evidence for robust circulation. The expectation for BD spectra at temperatures at or below the L/T transition is to favor higher relative abundances of CH₄ and NH₃ and lower abundances for CO and N₂. However, numerous BD atmospheric observations dating back to Gliese 229B have detected above equilibrium values of CO (Noll et al. 1997; Oppenheimer et al. 1998) and a dearth of NH₃ (Saumon et al. 2000, 2006).

This disequilibrium was expected by atmospheric modelers based on observations of Jupiter and Saturn (Fegley & Lodders 1996). Chemical disequilibrium might indicate atmospheric circulation with convection or advection transporting CO into the higher stratosphere while sequestering NH₃ within the BD and gas giant interiors (Noll et al. 1997; Griffith & Yelle 1999; Saumon et al. 2000; Leggett et al. 2007).

Applying Doppler imaging to CO, CH₄, and NH₃ spectral lines could potentially create BD and gas giant exoplanet surface maps detailing regions of increased convective circulation as well as other types of vertical transport observed in solar system stratospheres. Connecting these maps to dynamical models would allow for a greater understanding of ultracool objects' atmospheric structure and dynamics.